########################################################
## Interpretation of results from our PML estimators###
## To be run after Mlogit_NAPPHeterofformula1v2.R. To get Marginal Effects and every equilibria
## In 1. We compute ME for endogenous and exogenous effects
## In 2. We compute standard errors for marginal effects of the endogeneous effect via the Krinsky and Robb method
## In 3. we use spectral methods to compute every possible equilibria of our Nonlinear System
###################


rm(list=ls(all=TRUE))
set.seed(12345)
root <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/" #root <- "C:/HPC/"
#root <- "C:/HPC/"
setwd(paste(root,"Results/ReStat/",sep=""))
#codes <- "C:/HPC/CODE/"
codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/C_ANALYSIS/" #codes <- "C:/HPC/CODE/"
source(paste(codes,"Mlogit_NAPPHeteroFunctionsv2.R",sep=""))
foo(c("BB","mvmeta")) #mvmeta for cbindlist
# We have to change the following line depending on whether we are using the baseline exercise or a robustness check.
# notice that the ending of the dataset is changed manually
wexernei <- "neww50list" #"neww100list" "newnmovers" "newage1530" "newnatleast1"
AM <- 0 #if using Aguirregabiria and Mira==1, 0 if using Kasahara and Shimotsu
if( AM==1 ){  rootsave <- paste(root,"Results/ReStat/",sep="") } else   {rootsave <- paste(root,"Results/ReStat/KS/",sep="")}
sink(paste(rootsave,"NAPPEstimationHetewfeInterpret",wexernei,".txt",sep=""))

#load original dataset
load(file=paste(root,"NAPP/ReStat/pooleddataNAPP",wexernei,".RData",sep="" )) # pooleddata, Data per parish in long format (sw) pooled together
load(file=paste(root,"NAPP/ReStat/pooleddataNAPPhomo",wexernei,".RData",sep="" )) # pooleddatahomo, Data per parish in long format (sw) pooled together for hmogenous
load(file=paste(root,"NAPP/ReStat/pooleddatalistNAPP",wexernei,".RData",sep="" )) #pooleddatalist, Data per parish in wide format (sw.1,sw.2,...) pooled together
load(file=paste(root,"NAPP/ReStat/pooleddatalistNAPPhomo",wexernei,".RData",sep="" )) #pooleddatalisthomo, Data per parish in wide format (sw.1,sw.2,...) pooled together
load(file=paste(root,"NAPP/ReStat/datalistNAPP",wexernei,".RData",sep="" )) #datalist, Lists of data per parish in wide format (sw.1,sw.2,...)
load(file=paste(root,"NAPP/ReStat/coordslistNAPP",wexernei,".RData",sep="" )) #coordslist, Lists of coordinates per parish
load(file=paste(root,"NAPP/ReStat/dneighboursNAPP",wexernei,".RData",sep="")) #dneighbours, Lists of neigbours identities per parish
load(file=paste(root,"NAPP/ReStat/wlistNAPP",wexernei,".RData",sep="")) #wlist, Sparse Lists of neighbours per parish
data <- split(pooleddata, list(pooleddata$IDSocial)) #data, Lists of data per parish in long format (sw)
#Load last PML results
load(paste(rootsave,"JePMLmtwfinal",wexernei,".RData",sep=""))
load(paste(rootsave,"PMLmtwfinal",wexernei,".RData",sep=""))

###
#1#
###
altern <- sort(unique(pooleddata$alt))
#1.1To compute ME from the endogenous effect. dP_{iy}/dw_is_y0= J_y P_{iy} (1-P_{iy}) if y0==y; -J_y0 P_{iy0} P_{iy}
if (ncol(ml.m.w$probabilities)<length(c(min(altern):max(altern)))) {
prob <- cbind(apply(ml.m.w$probabilities,1,function(x) 1-sum(x)),ml.m.w$probabilities)
} else prob <- ml.m.w$probabilities 
colnames(prob) <- c(min(altern):max(altern))
MEb<-apply(prob,2,function(x) x * (1-x))
ME <- t(t(MEb) * summary(ml.m.w)$coeff[paste("sw:",c(min(altern):max(altern)),sep="")])
print("Marginal Effects:")
print(apply(ME,2,mean))

#To give some sense, compute sd of sw.x (in equilibrium) and postmultiply such ME
for (alter in c(min(altern):max(altern))) {
#if (alter==0) sdswe <- c(sd(ml.m.w$data[ml.m.w$data$alt==alter,]$sw)) else  sdswe <- c(sdswe,sd(ml.m.w$data[ml.m.w$data$alt==alter,]$sw))
  if (alter==0) sdswe <- c(sd(ml.m.w$model[ml.m.w$model$alt==alter,]$sw)) else  sdswe <- c(sdswe,sd(ml.m.w$model[ml.m.w$model$alt==alter,]$sw))
if (alter==0) sdsw <- c(sd(pooleddata[pooleddata$alt==alter,]$sw)) else  sdsw <- c(sdsw,sd(pooleddata[pooleddata$alt==alter,]$sw))

}
print("std dev sw estimated:")
print(sdswe)
print("Marginal Effects * std dev sw estimated:")
print(apply(ME,2,mean)*sdswe)

print("std dev sw:")
print(sdsw)
print("Marginal Effects * std dev sw:")
print(apply(ME,2,mean)*sdsw)

#1.2To compute ME from the exogenous effect for continuous vars, see formula in our paper
Wlist <- lapply(wlist, as_dgRMatrix_listw) #converting neighbours list into a weighting matrix notice the W in upper case
Wlistprov <- lapply(Wlist,function(x) x * 1) #to convert it to a dgCMatrix

cvars <- c("AGE","nchild","nservant","married","stayerp") #cvars <- c("married","stayerp")
for (cvar in cvars) {
#cvar <- "AGE"
  cks <- summary(ml.m.w)$coeff[c(sapply(cvar, function(x) paste(x,":",c(min(altern):max(altern)),sep="")))]
  cks[is.na(cks)]<-0
  dks <- summary(ml.m.w)$coeff[c(sapply(cvar, function(x) paste(x,"w:",c(min(altern):max(altern)),sep="")))]
  dks[is.na(dks)]<-0
  Js <- summary(ml.m.w)$coeff[paste("sw:",c(min(altern):max(altern)),sep="")]
  ##1.2.1. direct effect
  cksprob <- rowSums(prob %*% diag(cks))
  for (cl in c(1:length(cks))) {
    MEdirect <-apply(prob,2,function(x) x * (cks[cl]- cksprob))  
  }
  print(paste("direct:",cvar))
  #print(apply(MEdirect,2,mean)*sd(pooleddata[[cvar]]))
  print(apply(MEdirect,2,mean))
  
  ##1.2.2. indirect through y
  dksprob <- rowSums(prob %*% diag(dks))
  for (cl in c(1:length(dks))) {
  MEindirecty <-apply(prob,2,function(x) x * (dks[cl]- dksprob))  
  }
  MEindirecty <- cbind(pooleddatalist$IDSocial,MEindirecty)
  colnames(MEindirecty)[1]<-"IDSocial"
  MEindirecty <- as.data.frame(MEindirecty)
  MEindirectylist <- split(MEindirecty,list(MEindirecty$IDSocial))
  #MEindirectylist <- lapply(MEindirectylist,function(x) x[,c(2:6)]) #keeping only the probabilities
  MEindirectylist <- lapply(MEindirectylist,function(x) x[,c((2+min(altern)):(2+max(altern)))]) #keeping only the probabilities 
  multi <- function(x,y) {
    z <- x * y[[yy]] 
    return(z)
  }

  multimat <- function(x,y) {
    z <- diag(x %*% y)
    return(z)
  }
  WMEindirectylist <- vector(length=length(colnames(prob)),mode="list")
  ii <- 1
  for (yy in colnames(prob)) {
    prov <- mapply(multi,Wlist,MEindirectylist) 
    
    WMEindirectylist[[ii]] <- mapply(multimat,Wlistprov,prov)  
    WMEindirectylist[[ii]] <- unlist(WMEindirectylist[[ii]])
    ii <- ii+1
  }

  cbindlist <- function(list) {
    n <- length(list)
    res <- NULL
    for (i in seq(n)) res <- cbind(res, list[[i]])
    return(res)
  }

  WMEindirectylist <- cbindlist(WMEindirectylist)
  MEIndirecty <- ME * WMEindirectylist
  
  ##1.2.3. indirect through all y^0 \neq y
  MEIndirecty0prov <- t(t(prob * WMEindirectylist) * (-1)*(summary(ml.m.w)$coeff[paste("sw:",c(min(altern):max(altern)),sep="")]))
  MEIndirecty0 <- MEIndirecty0prov
  ii <- 1
  for (yy in colnames(prob)) {
    positions <- setdiff(c(1:length(colnames(prob))),ii)
    MEIndirecty0[,ii] <- prob[,ii]  * rowSums(MEIndirecty0prov[,positions])
    ii <- ii+1
  }


  sdvar <- sd(pooleddata[[cvar]]) 
  print(paste("std dev",cvar))
  print(sdvar)

  print(paste("Marginal direct Effect * std dev",cvar))
  print(apply(MEdirect,2,mean)*sdvar)
  print(paste("Marginal indirect y Effects * std dev",cvar))
  print(apply(MEIndirecty,2,mean)*sdvar)
  print(paste("Marginal indirect all y0 Effects * std dev",cvar))
  print(apply(MEIndirecty0,2,mean)*sdvar)
}
sink()

###
#2#
###
#################################################
#To compute standard errors for marginal effects of the endogeneous effect via the Krinsky and Robb method
# See Dowd-Greene-Norton2014_Computation of Standard Errors
#################################################
#See for some extended functions for mnlogit objects https://github.com/cran/mnlogit/blob/master/R/mnlogit.methods.R 

rm(list=ls(all=TRUE))
set.seed(12345)
root <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/" #root <- "C:/HPC/"
#root <- "C:/HPC/"
setwd(paste(root,"Results/ReStat/",sep=""))
#codes <- "C:/HPC/CODE/"
codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/C_ANALYSIS/"
#codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Codes/RCodes/EEACodes/" #codes <- "C:/HPC/CODE/"
source(paste(codes,"Mlogit_NAPPHeteroFunctionsv2.R",sep=""))
source(paste(codes,"mnlogit.methods.R",sep=""))
foo(c("micEcon", "mvtnorm")) #"micEcon" to use symMatrix, "mvtnorm" to use rmvnorm
foo(c("BB","mvmeta")) #mvmeta for cbindlist
wexernei <- "neww50list" #"neww100list" "newnmovers" "newage1530" "newnatleast1"
AM <- 0 #if using Aguirregabiria and Mira==1, 0 if using Kasahara and Shimotsu
if( AM==1 ){  
  rootsave <- paste(root,"Results/ReStat/",sep="") 
  } else   {rootsave <- paste(root,"Results/ReStat/KS/",sep="")}

sink(paste(rootsave,"AMEStdErrors",wexernei,".txt",sep=""))

S <- 1000 # Number of replications
load(paste(rootsave,"PMLmtwfinal",wexernei,".RData",sep="")) # to get estimates from the PML with parish fe
vcovmat <- vcov.mnlogit(ml.m.w) #to get variance covariance matrix
#to convert variance covariance matrix into a symmetric one (we know it is symmetric due to machine precision the matrix is never so)
vcovmat <- symMatrix(vcovmat[upper.tri(vcovmat, TRUE)], nrow=nrow(vcovmat), byrow=TRUE) 

#to get original data from estimation model
mlmwdata <- ml.m.w$model[ order(as.numeric(row.names(ml.m.w$model))),]
mlmwdatamlogit <- mlogit.data(mlmwdata,choice="choice",shape="long",alt.var="alt")

# to get estimated coefficientes
coeff <- summary(ml.m.w)$coeff
#to get S samples of the multivariate normal with mean coeff and variance covariance vcovmat to apply KR method
betas <- rmvnorm(n=S,mean = coeff,sigma = vcovmat) 
L <- length(c(min(mlmwdata$alt):max(mlmwdata$alt)))

AME <- matrix(rep(0,S*L),nrow=S) #Matrix to keep track of AME for each s.
AMEswe <- matrix(rep(0,S*L),nrow=S) #Matrix to keep track of AME*swe for each s.

ml.m.wKR <- ml.m.w

#To give some sense, compute sd of sw.x (in equilibrium) and postmultiply such ME
for (alter in c(min(mlmwdata$alt):max(mlmwdata$alt))) {
  if (alter==0) sdswe <- c(sd(ml.m.w$model[ml.m.w$model$alt==alter,]$sw)) else  sdswe <- c(sdswe,sd(ml.m.w$model[ml.m.w$model$alt==alter,]$sw))
}
print("std dev sw estimated:")
print(sdswe)

#2.1 To compute ME from the endogenous effect. dP_{iy}/dw_is_y0= J_y P_{iy} (1-P_{iy}) if y0==y; -J_y0 P_{iy0} P_{iy}
s <- 1
while (s <= S) {
  print(s)
  ml.m.wKR$coefficients <- betas[s,]# Change coefficients beta_s
  prob <- predict(ml.m.wKR, newdata=mlmwdatamlogit) # predict probabilities
  MEb<-apply(prob,2,function(x) x * (1-x)) # compute marginal effects
  ME <- t(t(MEb) * summary(ml.m.wKR)$coeff[paste("sw:",c(min(mlmwdata$alt):max(mlmwdata$alt)),sep="")]) 
  AME[s,] <- apply(ME,2,mean,na.rm=TRUE)
  AMEswe[s,] <- AME[s,]*sdswe 
  s <- s+1
}

# Computing the SE via Krinsky and Robb method
AMEmean <- apply(AME,2,mean,na.rm=TRUE)
print("Mean AME KR method")
print(AMEmean)

AMEse <- apply(AME,2,sd,na.rm=TRUE)
print("std error of AME")
print(AMEse)

AMEswemean <- apply(AMEswe,2,mean,na.rm=TRUE)
print("Mean of AME*swe KR method")
print(AMEswemean)

AMEswese <- apply(AMEswe,2,sd,na.rm=TRUE)
print("std error of AME*swe")
print(AMEswese)

sink()

# ###
# #3#
# ###
# #To compute every possible equilibria per parish
rm(list=ls(all=TRUE))
set.seed(12345)
root <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/" #root <- "C:/HPC/"
#root <- "C:/HPC/"
setwd(paste(root,"Results/ReStat/",sep=""))
#codes <- "C:/HPC/CODE/"
codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/C_ANALYSIS/"
#codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Codes/RCodes/EEACodes/" #codes <- "C:/HPC/CODE/"
source(paste(codes,"Mlogit_NAPPHeteroFunctionsv2.R",sep=""))
wexernei <- "neww50list" #"neww100list" "newnmovers" "newage1530" "newnatleast1"
AM <- 0
if( AM==1 ){  rootsave <- paste(root,"Results/ReStat/",sep="") } else   {rootsave <- paste(root,"Results/ReStat/KS/",sep="")}

#load original dataset
load(file=paste(root,"NAPP/ReStat/pooleddataNAPP",wexernei,".RData",sep="" )) # pooleddata, Data per parish in long format (sw) pooled together
load(file=paste(root,"NAPP/ReStat/pooleddataNAPPhomo",wexernei,".RData",sep="" )) # pooleddatahomo, Data per parish in long format (sw) pooled together for hmogenous
load(file=paste(root,"NAPP/ReStat/pooleddatalistNAPP",wexernei,".RData",sep="" )) #pooleddatalist, Data per parish in wide format (sw.1,sw.2,...) pooled together
load(file=paste(root,"NAPP/ReStat/pooleddatalistNAPPhomo",wexernei,".RData",sep="" )) #pooleddatalisthomo, Data per parish in wide format (sw.1,sw.2,...) pooled together
load(file=paste(root,"NAPP/ReStat/datalistNAPP",wexernei,".RData",sep="" )) #datalist, Lists of data per parish in wide format (sw.1,sw.2,...)
load(file=paste(root,"NAPP/ReStat/coordslistNAPP",wexernei,".RData",sep="" )) #coordslist, Lists of coordinates per parish
load(file=paste(root,"NAPP/ReStat/dneighboursNAPP",wexernei,".RData",sep="")) #dneighbours, Lists of neigbours identities per parish
load(file=paste(root,"NAPP/ReStat/wlistNAPP",wexernei,".RData",sep="")) #wlist, Sparse Lists of neighbours per parish
data <- split(pooleddata, list(pooleddata$IDSocial)) #data, Lists of data per parish in long format (sw)
#Load last PML results
load(paste(rootsave,"JePMLmtwfinal",wexernei,".RData",sep=""))
load(paste(rootsave,"PMLmtwfinal",wexernei,".RData",sep=""))
secluster <- c(7.7274e-01,7.8952e-01,5.5157e-01,8.1746e-01,6.3272e-01,2.6553e+00) #clustered std errors of our main specification
sink(paste(rootsave,"Equilibria.txt",sep=""))
options(max.print=1000000)
ml.m.w0 <- ml.m.w
exer <- 1
while (exer <3){
  print("Multiple equilirbia")
  ml.m.w <- ml.m.w0
  if(exer ==1) { 
    print("Original coeff")
  } else {
      print("Original coeff + 1se")
    ml.m.w$coefficients[paste("sw:",c(min(ml.m.w$model$alt):max(ml.m.w$model$alt)),sep="")] <- ml.m.w$coefficients[paste("sw:",c(min(ml.m.w$model$alt):max(ml.m.w$model$alt)),sep="")]+secluster
    }
  
  # source(paste(codes,"Mlogit_NAPPHeteroFunctions.R",sep="")) #code to recover sw
  foo(c("BB"))
  source(paste(codes,"Mlogit_MontecarloFunctions.R",sep="")) #code to produce pparseFormula
  source(paste(codes,"roots.social.mnlogitv1.R",sep="")) #code to produce ProbMat
  #probMat<- roots.social.mnlogit(ml.m.w, newdata=data[[1]], probability=TRUE,wdata=wlist[[1]],p0=cbind(1-rowSums(ml.m.w$probabilities),ml.m.w$probabilities))
  # 
   genstart <- function(np,K) { #number of observation, number of alternatives, number of starting values
     p0 <-  matrix(runif(np*K),np,K)
     p0 <- c(matrix(p0/rowSums(p0), ncol = 1, byrow=FALSE))
     p0
   }
   genstartext <- function(np,K) { #number of observation, number of alternatives, number of starting values
     p01 <-  rep(1,np)
     p00 <- rep(0,(np*K-np))
     p0 <- matrix(c(p01,p00),np,K)
     p0 <- t(sapply(c(1:K),function(x) p0[,c(setdiff(c(1:K),c(1:x)),c(1:x))] ))
     p0
   }
   #Function to use inside BBsolve
   roots.social <- function(s) {
     n <- length(s)
     r <- s - roots.social.mnlogit(ml.m.w, newdata=data[[pp]], probability=TRUE, 
                                   wdata=wlist[[pp]],p0=s)
     r
   }
   
   #Doing iteration for parsihes
   P <- length(data) #total number of parishes
   nms <- 1000 #number of starting values
   #nms <- 2 #number of starting values
   equilibria <- vector(length=P,mode="list") #keep track of equilibria per group
   for (pp in c(1:P) ) {
     print(paste("orden data: ",pp,sep=""))
     print(paste("IDSocial: ",unique(data[[pp]]$IDSocial),sep=""))
     K <- 6 #length(unique(pooleddata$alt))
     np <- nrow(data[[pp]])/K ## number of members in parish 
     #initial probabilities, single roots
     #p0 <-  matrix(runif(np*K),np,K)
     #p0 <- c(matrix(p0/rowSums(p0), ncol = 1, byrow=FALSE))
     #root<- BBsolve(par = p0, fn = roots.social)$par
     
     #initial probabilities, multiple roots
     p0 <-  t(sapply(c(1:nms),function(x) genstart(np,(x-x)+K))) #each row has a set of starting values
     p0 <- rbind(p0,genstartext(np,K)) #attaching together the extreme starting values, where in each row inds follow each occupation with prob 1
     ans <- multiStart(par = p0, fn = roots.social, action = "solve")
     sum(ans$conv)  
     pmat <- ans$par[ans$conv, ] 
     ord1 <- order(pmat[, 1])
     ans <- round(pmat[ord1, ], 4)
     equilibria[[pp]] <-ans[!duplicated(ans), ] 
     print(paste("Number of Equilibria: ",sum(!duplicated(ans)),sep=""))
   }
exer <- exer+1
}
 sink()
